Integrable two-channel + ipy-wave superfluid model 
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We present a new two-channel integrable model describing a system of spinless fermions interacting 
through a p-wave Feshbach resonance. Unlike the BCS-BEC crossover of the s-wave case, the p- 
wave model has a third order quantum phase transition. The critical point coincides with the 
deconfinement of a single molecule within a BEC of bound dipolar molecules. The exact many- 
body wavefunction provides a unique perspective of the quantum critical region suggesting that the 
size of the condensate wavefunction, that diverges logarithmically with the chemical potential, could 
be used as an experimental indicator of the phase transition. 
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In recent years p-wave paired superfluids have at- 
tracted a lot of attention, in part due to their exotic 
properties Of particular interest is the chiral two- 
dimensional {2D) Px + ipy superfluid of spinless fermions, 
that supports a topological phase with zero energy Ma- 
jorana modes 0]- The latter are theorized to serve as 
a basic element for a topological quantum computer Q. 
That exotic superfluid state might be realized in the Al 
phase of ^He [J], in the layered perovskite oxide Sr2Ru04 
[5|, and in the Pfafhan quantum Hall state at v = 5/2 
filling Q. Most promising is its realization in a cold 
gas of fermionic atoms in a single hyperfine state. In- 
deed, p-wave Feshbach resonances have been observed 
and studied in ^Li and Q with the potential to ma- 
nipulate the system from the weak (BCS) to the strong 
(BEC) pairing regime. However, these gases revealed to 
be unstable due to atom-molecule relaxation processes in 
which the molecule decays to a deep bound state while 
the atom escapes with excess energy Q. Other atomic 
and molecular gases are now in consideration, as well as 
different mechanisms to suppress relaxation. 

From a theoretical standpoint, despite great efforts to 
describe these systems, a complete understanding of the 
BCS-BEC transition and the corresponding phase dia- 
gram is still missing. Recently, by means of an exactly 
solvable p^ + ipy pairing model 0, [l^l it was shown that 
the quantum phase transition (QPT) taking place from 
weak pairing to strong pairing can be understood as the 
deconfinement of bound Cooper pairs [lH. In this Let- 
ter we introduce an exactly solvable two-channel Px + ipy 
pairing model with a Feshbach resonance. Our model 
is exactly solvable in arbitrary dimensions, although we 
will concentrate on its 2D realization. We propose a way 
to experimentally detect the so-called topological QPT 
[l^ from weak to strong pairing, by measuring a density- 
density correlation function. This together with the anal- 
ysis of the size of a Cooper pair in terms of the exact 
solution allows the characterization of the transition as 
one of a confinement-deconfinement type without Lan- 
dau order parameter. Moreover, the transition is shown 



to be third order in the Ehrenfest classification. Another 
way to theoretically detect that QPT is by analyzing the 
behavior of the quantum fidelity z of the ground state 
wavefunction. Interestingly, the second order derivative 
of ln|z| displays a logarithmic singularity at the transi- 
tion point confirming its third order character. 

Consider the 2D two-channel px + ipy-wave model 
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where cj^, creates a fermion in mode k 
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is a bosonic creation operator, and Hb = 6b^b + 
' h^Wbb. Our next goal is to show that this model is 
a particular realization of a family of exactly-solvable 
atom-molecule Hamiltonians of physical relevance in the 
context of cold atom physics. 

We start our demonstration by recalling the integrals 
of motion of the hyperbolic Richardson-Gaudin model 
which can be generically written as 
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where Sf, Sf, are the three generators of the SU{2)i 
algebra of mode i, i = 0, • • • , L, with spin representation 
Si such that (S'f) = Si{si~\-1). Therefore, the operators i?i 
contain L free parameters rji plus the strength of the 
quadratic term 7. The integrals of motion ([2]) commute 
with the z component of the total spin component 5*^ = 

We next single out the copy i = and consider its 
large spin sq limit; eventually we are interested in the 
limit So ~^ 00. The corresponding SU{2)q generators are 
bosonized by means of the Holstein-Primakoff mapping 



S'o — 



So, 



S+ = b^ V2sn ~ Wh, So = (5'o+)^ 
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In the spin-boson representation the conservation of 
becomes b^b + J2i=i = M — Lc/2, where Lc = 

Inserting the boson representation into the integrals of 
motion and expanding them in terms of I/sq, we arrive 
to the complete set of integrals of motion describing a 
spin-boson model 
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where we took advantage of the freedom to select the 
values of 7 and 770: 7 = — 1- -r^ and 770 = 2g^so, so 
that finite integrals of motion result in the limit sq — > 00 . 
The detuning parameter S in Hh is given hy S = {Lc — 
2{M — 1) — k)(7^/2, with k and free parameters. 
The corresponding eigenvalues [l^ in this limit are 
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FIG. 1: Pairon distribution for x = —0.1 (triangles), x = 0.37 
(solid circles) and x — 1.7 (squares) in a system of M = 10 
pairons, Lc = 40 [p — 1/4), A = g^Lc = 1/2, and uj = 1. 
The open circles represent the five lowest level parameters 
j/fe. The upper inset shows the bosonic, po, and fermionic, 
p — po, densities for the exact solution (solid line) and the 
BCS approximation (dashed line). The lower inset displays 
the exact momentum distribution, , for the three cases and 
the BCS approximation in solid lines. 
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where the set of pair energies (pairons) Eg represents a 
particular solution of the Richardson equations 
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with Qo = (M - 1) - Lc/2 + i6/2g^). The exact eigen- 
states in turn are given by 
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where |0, z.) represents the boson vacuum tensor and 
the state of seniority ly, such that Sf |^) = — ji^) and 
S- \v) = for all i. 

This completes the derivation of the integrable spin- 
boson model defined by the integrals of motion ([3]) , whose 
eigenvalues and common set of eigenvectors arc given by 
(|4]) and ([6|) respectively, expressed in terms of the pair 
energies Ea solutions of the Richardson equations ([5]). 

The connection between the spin-boson model above 
and the + ipy pairing model of Eq. ([1]) is realized 
by the pair representation of the SU{2) algebra = 

^{hk + n-k-1), = {S~y = {ka: + iky)clc'_i^/ \k\, 
where now we consider the mode index to represent mo- 
mentum k = (kx, ky). Now M represents the total num- 
ber of fermionic pairs and bosons and Lc is the maxi- 
mum possible number of fermionic pairs. Inserting this 



representation in the integral of motion TZq and defin- 
ing r/k ~ we obtain the Hamiltonian ([1]), H = TZq, 
thus showing that it is exactly solvable. Indeed, any lin- 
ear combination of the integrals of motion defines an ex- 
actly solvable Hamiltonian. In particular, the Hamilto- 
nian studied by Links et al. [l5|, can be obtained from 
the linear combination H = 2'^^eiTZi with 2ei = l/rji. 
If efe = |fcp/2, with ^ /c ^ fccut, then the parameters 
rjk entering in the Richardson's equations ([5]) are defined 
in the interval l/j/ccutP — l/(2a;) ^ 77^ ^ 00. Their 
exactly solvable model does not display a QPT between 
the weak and strong pairing phases, because, as stated in 
Ref. llj , for the Px + ipy model to display a non-analytic 
behavior in the continuum limit it is required that one of 
the parameters rjk vanishes for a given mode (e.g., k ~ 
mode). 

Coming back to the Eqs. in complete analogy with 
Ref. |ll| . we can recognize two special cases. One in 
which all the pairons Ea converge to zero, corresponding 
to 2Q0 + I = AI, and a second case where the first pairon 
converges to zero, i.e., when 2Qo -1-1 = 0. Defining x = 
S/{2g'^Lc), these two cases correspond to 2a:MR = 1 — 
{M - l)/Lc and 2xcr = 1 - (27\/ - l)/Lc. The first case 
defines the so-called Moore-Read line. We will show that 
the second case signals an interesting third-order QPT. 

To get insight into the properties of the different phases 
we analyze the behavior of the pairons in a finite system. 
Figure [1] displays the pairon distribution for x < Xcr 
(strong pairing), Xcr < x < xmr and xmr < x (weak 
pairing on both sides of the Moore-Read line) for a sys- 
tem consisting of M — 10 pairs lying in a disk of radius 
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five units in an otherwise square lattice with Lc = 40 
[p = M/Lc = 1/4) and u = 1. For x = 1.7 the sys- 
tem has a fraction of 6 Cooper pairs (complex pairons) 
and 4 quasi-free pairs states (almost real and positive 
pairons). Crossing the Moore- Read line, at x = 0.37, the 
system is a mixture of 8 Cooper pairs and 2 bound pairs 
(real an negative pairons) . Well inside the strong pairing 
phase, at X = —0.1, all pairs are bound. Therefore, the 
weak pairing phase is characterized by a mixture of free 
fermions. Cooper pairs and bound molecules, while the 
strong pairing phase is a BEC of bound molecules. The 
upper inset displays the bosonic density, po = /Lc, 
while the lower inset shows the major rearrangement that 
takes place in the momentum distribution Uk close to 
fc = between the weak and the strong pairing phases. 
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FIG. 2: Second (inset) and third-order derivatives of the en- 
ergy density e in the thermodynamic limit (p = 1/4, A — 1/2, 
and to = 1 leading to Xcr = 1/4). 

Let us now consider the thermodynamic limit {L^ — > oo 
with A — g^Lc finite) of this two-channel p-wave model. 
Following the same procedure as in [llj we obtain a pair 
of couple BCS-like equations for the unknowns pa and 
chemical potential p 
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uj - p + Xpa + f 
Xpa - p + \p\ 

uj - p + Xpo + f 
Xpo- p + \p\ 



(7) 
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where / = ^J 2Xujpo + {uj — p)^ . 

The ground state energy density e = E/Lc is 

. , s , > f 1 , A ^Apo p+ \p\ 
e = p(uj + p) + Xpox ~ [^^ 2 J ~ 2 — ■ 

If we identify Apo = ^Ics ^^"^ ^ ~ I/S^bcs, the energy 
density e coincides with the energy density for the one- 
channel P:c + ipy model . Analogously, the right hand 
sides of Eqs. (I7][S|) coincide with the corresponding ones 
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FIG. 3: Quantum phase diagram in terms of p and x, for a 
fixed value of A. Phase boundaries are insensitive to A. 



in However, their left hand sides arc different due 
to the existence of two independent parameters x and A, 
or equivalently 6 and g. Potential non-analyticities at 
p = can be attributed to the presence of absolute value 
of p terms. Indeed, from Fig. [5] that shows the second- 
and third-order derivatives of the energy density with 
respect to the control parameter x, we conclude that the 
p-wave atom-molecule model displays a third-order QPT 
at the critical value Xcr = — p, coinciding with the 
limit in which all pairons are real and negative except 
one that converges to zero, and p = 0. The Moore-Read 
line corresponds to xmr = (1 ~ p)/2 and p = Apo/2, 
leading to e = 0, coinciding with the limit in which all 
pairons converge to zero. The resulting quantum phase 
diagram depicted in Fig. [3J depends on the density p, and 
on the control parameters A and x. However, the phase 
boundaries are independent of A. Thus, our two-channel 
p-wave model extends the one-channel model to: a) x < 
due to the possibility of having negative detunings S, b) 
p > 1 due to the mixture with a bosonic system, and c) 
include an extra control parameter A. 

Since there is no Landau order parameter characteriz- 
ing this third-order transition, one would like to devise 
a way to experimentally detect it. In Ref. the root 
mean square of the condensate wave function was pro- 
posed as a possible indicator of that QPT 
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where a and b are the 



V(|fcP-a)(|fc|2-6)' 

roots of the polynomial = {ek — pY + '^^Po^k = 
(efc — a/2)(efc — 6/2). The integral ^ can be obtained 
analytically. As can be seen in the right inset of Fig. 
m Trms diverges logarithmically at the QPT point as 
rLs -> - 9^'" w/i^t lr-^- Making contact with the ex- 

rms 2Apn ln(l+cj/2Apo) ^ 

act wave function ([6]) we can immediately associate this 
divergence to the confinement of the last Cooper pair 
or, coming from strong pairing, to the deconfinement 
of a single bound pair within the BEC molecular wave 
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function. Away from the QPT point, r^ms diverges again 
in the extreme weak couphng hmit, a; — > (X), in which 
all pairs are deconfined. The condensate wave function 
can be related to the density-density correlation function 
which in coordinate space can be expressed as 

\^ir - r')\' - ((n, - - {fir'))) + |i^(r - r')|^ 

where (f>{r — r') is the condensate wave function in coordi- 



nate space, Ur 



{firfi,.') is the density-density cor- 



relation function, and F{r — r') is the Fourier transform 
of the momentum density {fikj- In trapped cold atomic 
gases, the condensate wave function can be obtained from 
measurements of the density-density correlations using 
quantum noise intcrferomctry and the momentum dis- 
tribution from time-of-flight measurements after opening 
the trap. Therefore, the root mean square size of the 
condensate wave function constitutes a unique indicator 
of the QPT which can be experimentally accessed. 
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FIG. 4: Xf/Lc (see text) as a function of x, p — 1/4, A — 
1/2, oj = 1 for Ax = 10"^ (dot-dashed) , 5 x 10"" (solid), 
and 1Q~^ (dashed). At x^r = 1/4, as Ax — >■ 0, it diverges 
logarithmically. At x = Xcr ± Aa::, and finite Ax — 5 x 10~", 
the second-order derivative of x develops another logarithmic 
divergence in terms of fi. The root mean square of the Cooper 
pair size, rrms is shown in the right inset. 

The quantum fidelity z[x,Ax] = (*(.t - Ax}j^{x + 
Ax)) can also be used as an indicator of QPTs [16|. The 
quantity Xf = — 2 In |2;[x. Ax] j/Aa:^, which in the limit 
Ax — > is the so-called fidelity susceptibilty, has an 
essentially different behavior depending on the ground- 
state overlaps considered. If the overlap is taken between 
two states belonging to the same phase (ja; — Xcr] > 
Aa;), Xfl^c tends to a value which remains finite in 
the limit Aa; — ^ 0, whereas if the overlap is between 
two states in different phases (|a; — a;cr| < Ax) a In Aa; 
dependence appears. More precisely, {xf /Lc) \x^, ~ 
— {dfj./dx)^\n{Ax)/{ujXpo)\^ ■ Clearly, when Aa; — > 
the fidelity susceptibility develops a logarithmic diver- 
gence as a function of (x — Xcr) similar to that obtained 
in the one-channel px + ipy model [lol | . 

However, even for Aa; finite, Xfl^c shows a non- 
analytic behavior. When one of the ground-states in 
2;[a;, Aa;] is taken close to the transition point (a; w 



Xcr ± Ax), Xf/Lc can be written as {Xf / Lc)\^r^^^^j^^^ = 
-^|ln|/izp|/(2wApoAx2) + 0(^5},), with ^± = /i(x±Ax). 
At X = Xcr ± Ax, where /izp = /i(xcr) = 0, the previ- 
ous expression is continuous, but its second-order deriva- 
tive diverges logarithmically, as shown in Fig. 21 These 
results should be contrasted with the results obtained 



for an Ising chain of length L [l7|, where {xi / L)\^ ~ 
l/|Ax|, i.e. power law, and a logarithmic divergence at 
X = Xcr±Ax appears in the first-order derivative oixi I L. 

In conclusion, we presented a new exactly solvable 
spin-boson model which has as particular realization the 
two-channel px + ipy-wave superfluid. We showed that 
the model has a third order QPT that can be accessed 
experimentally by measuring the density-density corre- 
lation function. The analysis in terms of fidelity pro- 
vides further evidence of the non-Landau character of 
the QPT, and its logarithmic singularity indicates that 
it is of a confinement-deconfinement type, as suggested 
by the exact wave function. 
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